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1  Introduction 

1.1  Impact  and  Influence  of  the  j>- version 

Research  and  development  of  the  p-version  of  the  finite  element  method  at  Washington 
University  have  been  supported  by  AFOSR  since  1977. 

The  p-version  of  the  finite  element  method,  together  with  the  h—p  version,  have  now 
been  as  widely  accepted  as  the  classical  h-version,  as  eflficient  computational  approaches 
to  the  solution  of  a  wide  variety  of  engineering  and  scientific  problems.  One  current  com¬ 
puter  implementation  of  the  p-version  is  the  code  MSC/PROBE,  which  is  marketed  by 
the  MacNeal-Schwendler  Corporation,  p-version  capabilities  have  been  incorporated  into 
MSC/NASTRAN  and  are  now  available  to  the  general  user.  Thus,  the  p-version  which  until 
recently  has  only  been  licensed  for  use  by  a  large  number  of  aerospace  companies  and  govern¬ 
ment  laboratories  on  an  individual  basis,  is  now  accessible  to  ^  users  of  MSC/NASTRAN. 

PROBE  contains  many  unique  features  including  fracture  mechanics  extraction  proce¬ 
dures  which  were  developed  specifically  for  the  p-version  of  the  fimite  element  method.  These 
procedures,  along  with  explicit  quality  control  features,  provide  an  excellent  tool  for  perform¬ 
ing  fracture  mechanics  analyses  and  verifying  the  accuracy  of  the  results  [1]. 

Several  other  computer  codes  which  implement  the  p-version  are  now  available  for  engi¬ 
neering  use.  These  codes  include: 

•  ANSYS,  marketed  by  Swanson  Analysis,  Houston,  PA. 

•  IDEAS,  marketed  by  Structural  Dynamics  Research  Corp.,  Cincinnati,  OH. 

•  COSMOS/M,  marketed  by  the  Structural  Research  and  Analysis  Corporation,  Santa 
Monica,  CA. 

•  Mechanica,  marketed  by  Rasna  Corporation,  San  Jose,  CA. 

•  NISA/P-Adapt,  marketed  by  Engineering  Mechanics  Research  Corporation,  Troy,  MI. 

(for  some  details  on  these  codes,  see  Machine  Design,  July  25,  1991,  pp.  73-77,  and  Me¬ 
chanical  Engineering,  September,  1991,  pp.  83-84). 

Most  important,  perhaps,  is  a  new  software  structure,  PEGASYS,  marketed  by  ESRD 
(Engineering  Software  Research  and  Development). 

PEGASYS:  A  Software  Infrastructure  for  R  &  D 

PEGASYS  is  an  advanced  software  infrastructure,  designed  with  two  objectives  in  mind: 
To  support  research  and  to  provide  a  framework  for  technology  deployment  in  numerical 
simulation. 

Essentially,  PEGASYS  is  an  open  system  of  software  modules  with  well-documented  in¬ 
terface  specifications.  PEGASYS  currently  consists  of  approximately  250,000  lines  of  code, 
written  in  the  FORTRAN  and  C  languages.  PEGASYS  is  maintained  in  the  popular  work¬ 
stations,  such  as  Hewlett  Packard,  Silicon  Graphics,  SPARCstation,  IBM  RS6000  and  VAXs- 
tation.  Using  the  software  modules  of  PEGASYS,  it  is  possible  for  researchers  and  developers 
to  assemble  simulation  software  products  designed  to  meet  specific  engineering  and  scientific 
objectives. 
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A  number  of  finite  element  procedures  concerned  with  element  formulation,  mapping, 
assembly,  extraction  procedures,  error  estimation  and  adaptivity  have  been  completed.  The 
currently  active  R  &  D  projects  utilizing  PEGASYS  include  the  following: 

-  Development  of  hierarchic  models  for  structural  plates  and  shells; 

-  Development  of  hierarchic  models  for  laminated  composites; 

-  Superconvergent  methods  for  the  computation  of  stress  intensity  factors,  including 
generahzed  stress  intensity  factors  for  multi-material  interface  problems; 

-  Investigation  of  crack  propagation  phenomena; 

-  Mechanics  of  filament-wound  composites; 

-  Advanced  solution  methods  for  elastic-plastic  problems; 

-  Development  of  models  for  structural  connections; 

-  Error  estimation  and  quality  control  procedures; 

-  Implementation  of  advanced  finite  element  software  on  parallel  computers; 

-  Design  sensitivity  analysis. 

Some  current  R  &  D  users  of  PEGASYS  include:  Washington  University  in  St.  Louis, 
The  University  of  Maryland  at  College  Park,  The  University  of  Iowa,  McDonnell  Douglas 
Aerospace  Co.,  Ford  Motor  Co.,  Kelly  AFB,  The  University  of  Technology  in  Espoo,  Finland. 
PEGASYS  is  being  used  by  NASA  LBJ  Space  Center  in  Houston,  and  is  being  evaluated  for 
use  at  Lockheed  Fort  Worth,  McDonnell  Aircraft  Company  and  Ford  Motor  Company. 

Thus,  earher  research  into  the  p-version  which  was  supported  by  AFOSR,  has  rapidly 
been  transferred  into  modem  technology. 

1.2  Hierarchic  Modeling  and  Analysis  of  Structural  Connections 

In  the  last  few  years  there  has  been  a  great  deal  of  progress  in  the  development  of  methods 
for  controlhng  the  errors  of  discretization  in  finite  element  analysis.  It  has  been  estabhshed 
that,  with  properly  designed  finite  element  meshes  and  p-extensions,  it  is  possible  to  achieve 
exponential  convergence  rates  for  very  large  classes  of  problems  and,  with  superconvergent 
extraction  procedures,  engineering  data  can  be  computed  from  finite  element  solutions  with 
accuracies  comparable  to  the  accuracy  of  the  strain  energy.  The  stopping  criterion  is  that 
the  data  of  interest  must  he  substantially  independent  of  the  discretization. 

Control  of  errors  of  idealization  is  based  on  the  idea  that  any  particular  model,  for 
example  a  model  based  on  the  linear  theory  of  elasticity,  is  embedded  in  a  set  of  more 
general  models,  for  example,  ones  that  account  for  geometric  and/or  material  nonlinearities. 
Which  model  is  appropriate  for  a  given  task,  depends  on  the  goals  of  computation  and  the 
required  accuracy.  In  general  one  starts  with  a  simple  linear  model.  Once  a  reasonably 
accurate  approximate  solution  is  available,  then  it  is  possible  to  judge  whether  the  solution 
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violates  the  assumptions  (i.e.,  restrictions)  incorporated  in  the  model.  The  key  question 
is,  whether  removal  of  those  restrictions  would  have  a  significant  effect  on  the  conclusions 
drawn  for  the  model.  If  that  is  the  case  then  the  model  has  to  be  “upgraded”  and  a 
new  solution  obtained  which,  once  again,  has  to  be  examined  for  consistency  with  respect 
to  the  restrictions  incorporated  in  the  model.  The  stopping  criterion  is  that  the  data  of 
interest  must  be  substantially  independent  of  any  restrictions  imposed  on  the  model.  The 
construction  of  hierarchic  models  for  structural  plates  and  shells  is  described  in  reference  [2]. 
An  investigation  of  hierarchic  modelling  techniques  applied  to  fastened  structural  connections 
is  presented  in  reference  [3]. 

Fasteners  are  complicated  structural  assemblies,  full  stress  analysis  of  which  would  require 
consideration  of  three-dimensional  problems  involving  friction,  contact,  nonlinear  material 
properties,  and  other  factors,  such  as  mode  of  installation,  etc.  A  number  of  idealizing  de¬ 
cisions  are  necessary  in  order  to  make  the  problem  tractable  by  numerical  methods.  This 
generally  means  that  the  problem  has  to  be  reduced  to  a  two-dimensional  one  and  nonhnear- 
ities  must  either  be  ignored,  or  only  mild  forms  of  nonUnearities  considered.  The  objectives 
axe  to  compute  (a)  the  total  force  acting  on  fastener  groups;  (b)  the  forces  acting  on  individ¬ 
ual  fasteners;  (c)  critical  combinations  of  stresses  or  strains  in  the  neighborhood  of  the  most 
heavily  loaded  fasteners;  (d)  the  stress  intensity  factors  for  cracks  in  the  vicinity  of  fasteners. 

There  is  renewed  interest  in  this  problem,  particularly  in  the  aerospace  industry,  for  two 
reasons;  First,  concern  over  the  safety  of  aging  aircraft  requires  re-evaluation  of  structural 
connections  that  remain  in  service  much  longer  than  the  original  design  life.  Second,  previ¬ 
ously  established  design  procedures  are  based  on  assumptions  of  substantial  ductihty  in  the 
fasteners  and  the  connected  parts.  These  assumptions  do  not  hold,  in  general,  for  composite 
materials.  For  this  reason,  determination  of  the  distribution  of  forces  in  fastener  groups 
subjected  to  critical  loads  and  determination  of  the  stress  distributions  in  the  neighborhood 
of  fasteners,  are  important.  Given  the  very  low  ductihty  of  composite  materials,  use  of  the 
linear  theory  of  elasticity  is  sufficient  in  a  very  large  number  of  practical  problems. 

There  may  be  considerable  interest  on  the  part  of  the  Air  Force  in  analyzing  structural 
connections  because  there  is  a  potential  for  substantial  improvements  in  the  rehabihty  and 
performance  of  mathematical  models  for  structural  coimections  and  their  repairs.  Each 
repair  being  a  separate  design  problem,  improvements  in  this  area  are  of  obvious  importance 
to  the  Air  Force. 

The  representation  of  fasteners  in  current  modeling  practice  is  a  source  of  serious  errors: 
In  current  modehng  practice  fasteners  are  usually  represented  by  multipoint  constraints, 
that  is  nodes  of  the  finite  element  mesh  are  positioned  at  fastener  locations  and  the  nodes 
are  connected  by  rigid  or  fiexible  bars.  This  modeling  practice  is  conceptually  wrong  and  the 
computed  forces  in  the  fasteners  are  entirely  discretization-dependent  [3,  4,  5,  6,  7,  8]. 

1.3  References 

[1]  “Rehability  and  Quahty  Control  in  Fracture  Mechanics  Analysis  using  the  p- version  of 
the  Finite  Element  Program”,  MSC/Probe,  Proceedings  of  the  1990  USAF  Structural 
Integrity  Conference,  August  1991. 

[2]  “Hierarchic  Plate  and  Shell  Models  Based  on  p-Extension” ,  by  B.A.  Szabo  and  G.J. 
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2  Summary  of  Research  Accomplishments 

Need  for  iterative  solvers  and  parallelization  in  the  p-version. 

Direct  methods  axe  fuUy  satisfactory  for  the  solution  of  the  assembled  system  of  linear 
equations  in  the  p-version  for  small  size  problems  or  large  sparse  systems.  But  large  size 
problems  (with  more  than  30,000  degrees  of  freedom)  which  are  today  common  in  finite  ele¬ 
ment  analysis,  are  inefficient  both  in  terms  of  CPU  time  and  storage  [1].  Iterative  methods 
are  a  much  better  choice  for  problems  with  a  large  number  of  elements,  and  also  do  not  suffer 
from  fill-in.  Iterative  methods  can  be  implemented  on  parallel  computers  and  achieve  com¬ 
putational  load  balance,  that  is,  the  computational  load  can  be  balanced  among  processors. 
Classical  iterative  methods:  Jacobi,  Gauss-Seidel,  SOR,  SSOR,  Chebyshev  methods,  have 
been  studied  thoroughly  and  solve  well-conditioned  systems  efiiciently.  But  these  methods 
converge  very  slowly  for  ill-conditioned  systems. 

Also,  they  do  not  exploit  the  special  structure  of  the  global  stiffness  matrix.  Iterative 
methods,  specially  designed  for  the  p-version,  can  lead  to  dramatic  increases  in  the  rate 
of  convergence,  to  implementation  on  parallel  computers  and,  hence,  to  impressive  savings 
in  cost.  These  new  methods  will  be  used  for  complicated  problems  involving  hundreds 
of  thousands  of  degrees  of  freedom.  Such  problems  are  now  becoming  common  and  very 
reahstic. 

2.1  Parallel  Implementations  based  on  the  p-version 

An  iterative  method  based  on  the  textured  decomposition  method  has  been  developed  in 
order  to  solve  the  systems  of  linear  equations  arising  in  the  p-version  of  the  finite  element 
method.  The  iteration  was  used  to  implement  the  p-version  in  parallel  on  an  MIMD  computer 
NCUBE/six,  the  objectives  are  two-fold:  to  achieve  high  computational  efficiency  (that  is 
computational  load  should  be  balanced  among  the  processors)  and  simultaneously  to  achieve 
rapid  convergence. 

A  superelement,  consisting  of  four  adjacent  rectangular  finite  elements,  is  constructed  for 
two  dimensional  problems.  Based  on  the  structural  property  of  the  shape  functions,  each 
superelement  is  partitioned  into  three  blocks  in  two  different  ways,  and  a  two-leaf  textured 
decomposition  (TD)  is  used.  Computations  for  a  superelement  associated  with  each  leaf  are 
assigned  to  two  processors  and  are  performed  in  parallel.  A  new  preconditioner  is  introduced 
to  accelerate  convergence  in  a  preconditioned  textured  decomposition  (PTD).  A  special  local 
communication  strategy  is  used  to  avoid  global  assembly  and  global  communication. 

Three  model  problems:  A  Poisson  equation  on  a  rectangular  domain  with  a  smooth 
true  solution,  a  Laplace  equation  on  a  rectangular  domain  with  a  near  singular  solution, 
and  a  Poisson  equation  on  L-shaped  domain,  are  solved.  The  conjugate  gradient  method, 
the  textured  decomposition  method,  the  recursive  textured  decomposition  method,  both 
with  and  without  preconditioning;  and  the  classical  iterative  methods  (Jacobi,  Gauss-Seidel, 
SOR),  are  used  to  solve  the  three  model  problems.  Load  balance,  speedup  ratio,  and  spectral 
radii  of  the  various  iterations  are  studied.  The  test  results  indicate  that  recursive  PTD  with 
a  local  communication  strategy  gives  at  least  a  30%  improvement  in  computational  time 
over  the  other  methods. 
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This  method  has  been  presented  at  papers  at  the  SIAM  National  Meetings  and  will 
appear  in  the  SIAM  J.  on  Scientific  Computing  [2]. 

2.2  Multi-p  Methods 

A  natural  analogy  to  the  multigrid  method,  which  is  used  in  connection  with  the  finite  dif¬ 
ference  method  or  the  h-version  of  the  finite  element  method,  is  the  multi-p  method  which 
is  used  in  connection  with  the  p- version  of  the  finite  element  method  and  hieraxchical  shape 
functions.  Each  method,  multigrid  or  multi-p,  is  based  on  a  fundamental  iterative  scheme, 
e.g.,  Gauss-Seidel  or  SOR.  We  have  studied  multi-p  methods  based  on  general  linear  station¬ 
ary  iteration  schemes.  Various  V-cycle  algorithms  for  the  multi-p  methods  are  formulated 
including  standard  multi-p  V-cycle  (SMPV),  modified  multi-p  V-cycle  (MMPV)  and  vary¬ 
ing  multi-p  V-cycle  (VMPV).  Convergence  results  for  each  of  the  V-cycle  algorithms  have 
been  provided.  We  have  shown  that,  using  a  general  finear  stationary  iterative  scheme  as 
a  smoother,  the  standard  multi-  V-cycle  algorithm  has  a  linear  convergence  rate,  but  this 
rate  is  faster  than  that  of  the  smoother.  The  modified  multi-p  V-cycle  algorithm  also  has  a 
finear  convergence  rate,  but  again  this  rate  is  faster  than  the  rate  of  the  standard  multi-p 
V-cycle.  The  convergence  of  the  varying  multi-p  V-cycle  algorithm  is  a  consequence  of  the 
convergence  of  both  standard  and  modified  V-cycle  methods.  Numerical  experiments  on 
representative  problems  have  been  conducted,  and  the  numerical  results  agree  and  support 
our  theoretical  analysis  [3]. 

In  addition,  we  have  studied  nested  multi-p  methods.  An  error  estimate  has  been  de¬ 
rived  for  the  nested  multi-p  methods.  By  comparing  the  nested  multi-p  methods  with  the 
multi-p  V-cycle  methods,  we  found  that,  at  low  accuracy,  the  nested  multi-p  methods  are 
more  efficient,  but,  at  high  accuracy,  the  multi-p  V-cycle  methods.  This  leads  to  the  so- 
called  accelerated  multi-p  V-cycle  methods  which  are  a  combination  of  the  nested  multi-p 
methods  and  the  V-cycle  methods.  Convergence  of  the  accelerated  multi-p  V-cycle  methods 
is  proved.  Numerical  results  indicate  that  the  accelerated  multi-p  V-cycle  methods  are  80% 
more  efficient  than  the  underlying  iteration.  Some  of  the  numerical  results  are  shown  in  the 
figure  1. 
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(Preconditioned) 
Conjugate  Gradient 


Figure  1.  Comparison  of  the  multi-p  V-cycles. 

2.3  Multi-p  Preconditioning 

In  general,  preconditioned  conjugate  gradient  methods  are  regarded  as  very  promising  it¬ 
erative  methods  for  solving  linear  systems  of  equations.  In  particular,  conjugate  gradient 
methods  preconditioned  by  first  condensing  the  finite  elements  and  then  using  linear  elements 
to  construct  a  preconditioner  have  been  studied  by  some  researchers  [4]. 
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An  algebraic  theory  for  multi-p  methods  has  been  presented  and  analyzed.  Convergence 
and  symmetric  properties  are  proved  under  suitable  conditions.  It  is  then  shown  how  these 
multi-p  methods  can  be  used  as  preconditioners  for  the  conjugate  gradient  method  (CG). 
In  particular,  it  is  shown  that  given  any  preconditioner  Mp  to  CG,  a  multi-p  preconditioner 
Bp  based  on  Mp  can  be  constructed,  which  leads  to  a  smaller  condition  number  (and  hence 
faster  convergence).  When  Bp  is  applied  as  a  preconditioner  to  condensed  finite  elements, 
the  condition  number  is  shown  to  grow  slower  than  log^p),  the  best  currently  known 

result  for  2-D  problems  in  the  p- version  of  the  finite  element  analysis  [3]. 

Numerical  experiments  on  representative  problems  indicate  that  the  condition  numbers 
after  multi-p  preconditionings  axe,  in  fact,  independent  of  p.  The  numerical  results  also 
show  greater  efficiency  for  PCG  with  the  multi-p  preconditioners  in  terms  of  number  of 
iterations  and  CPU  time  when  compared  with  two  sophisticated  linear  equation  solvers:  (1) 
a  direct  frontal  solver  specially  designed  for  the  p-version  of  the  finite  element  analysis;  (2)  a 
highly  tuned  preconditioned  CG  code  in  ITPACK.  Preliminary  comparisons  of  the  number 
of  iterations  are  also  made  with  ROCKITS  [5],  a  new  commercial  code  used  for  the  p-version. 

The  methods  presented  are  intended  for  use  in  the  p-version  of  the  finite  element  analysis, 
but  are  general  in  nature  and  can  be  apphed  to  a  wide  variety  of  problems. 

The  following  three  dimensional  elastostatic  problem  on  a  brick  domain  illustrates  our 
results. 


Figure  2.  A  3-D  elastostatic  problem  on  a  brick  domain. 
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Figure  3.  Uniform  mesh  with  64  elements. 


A  3-D  elastostatic  problem  on  a  brick  domain-CPU  time  (sec.) 
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Figure  4.  Relationship  between  CPU  time  and  degree  p. 

These  results  have  been  reported  at  SIAM  National  Meetings  and  will  appear  in  SIAM 
J.  on  Scientific  Computing  [3,  6]. 

2.4  Conditioning  of  Global  Stiffness  Matrices 

In  order  to  properly  design  preconditioners  and  to  understand  the  rate  of  convergence  of 
iterative  Schemes,  it  is  important  to  know  the  condition  number  of  global  stiffness  matrices 
appearing  in  the  p-version. 

We  have  presented  a  theory  for  bounding  the  minimum  eigenvalues,  maximum  eigenval¬ 
ues  and  condition  numbers  of  stiffness  matrices  arising  from  the  p- version  of  finite  element 
analysis.  Both  lower  and  upper  bounds  are  derived  for  the  minimum  eigenvalues,  maximum 
eigenvalues  and  the  condition  numbers,  which  are  valid  for  stiffness  matrices  based  on  a  set  of 
general  basis  functions  that  can  be  used  in  the  p-version.  Although  the  theory  is  specifically 
developed  for  the  p-version  and  general  basis  functions,  it  is  shown  that  the  methodology 
used  can  be  apphed  easily  to  both  the  h- version  and  the  h  —  p  version  of  the  finite  element 
methods.  For  a  class  of  hierarchical  basis  functions  that  has  been  popularly  used  in  the  p- 
version,  explicit  bounds  are  derived  for  the  minimum  eigenvalues,  maximum  eigenvalues  and 
condition  numbers  of  stiffness  matrices.  Our  results  show  that  the  condition  numbers  of  the 
stiffness  matrices  grow  at  most  as  ,  where  d  is  the  number  of  dimensions.  This  growth 

is  quite  slow  compared  with  the  general  bounds,  which  provides  new  theoretical  support  for 
using  this  class  of  basis  functions.  Numerical  results  are  also  presented  which  indicate  that 
our  theoretical  bounds  are  quite  sharp. 
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2.5  The  Problem  of  Model  Selection 


This  phase  of  the  work  is  motivated  by  the  recognition  that  proper  model  selection  is  an 
essential  prerequisite  for  reliable  numerical  simulation  of  physical  systems.  Specifically,  work 
on  model  selection  is  focussed  on  two  areas:  One  is  the  development  of  proper  modelhng 
techniques  for  fastened  structural  connections,  the  other  is  model  selection  for  structural 
plates  and  shells  made  of  laminated  composites.  Both  areas  are  of  substantial  practical  im¬ 
portance:  In  computations  involving  durability  and  damage  tolerance,  refiable  and  accurate 
estimation  of  the  load  distribution  within  aeronautical  structures  is  essential.  Loads  are 
typically  transmitted  through  lugs,  actuators  and  fasteners.  Increasingly,  aircraft  compo¬ 
nents  are  fabricated  of  composite  materials.  Reliable  analytical  procedures  for  structural 
and  strength  analyses  of  these  important  structural  elements  are  not  currently  available. 

The  problem  of  fastened  structural  connections  is  approached  by  the  use  of  space  en¬ 
richment  techniques:  The  finite  element  space  is  enlarged  through  the  introduction  of  the 
fundamental  solution,  multiplied  by  a  cutoff  function,  such  that  the  nearly  singular  character 
of  the  fastener  interacting  with  the  plate  is  well  represented.  This  method  is  expected  to 
lead  to  convenient  and  accurate  representation  of  the  structural  action  of  large  numbers  of 
fasteners. 

In  the  case  of  laminated  plates  the  stress  distributions  near  boundaries,  discontinuities 
and  at  interfaces  are  generally  very  different  from  the  stress  distribution  in  the  interior  re¬ 
gions.  Boundary  layer  effects  are  normally  present,  and  the  problem  in  those  regions  is 
essentially  three-dimensional.  Hierarchic  models  for  laminated  plates  make  it  possible  to  ap¬ 
proximate  the  three-dimensional  problem  through  the  solution  of  two-dimensional  problems 
without  the  expense  of  a  fully  three-dimensional  analysis. 

The  proper  choice  of  model  depends  on  the  problem  description  and  the  data  of  interest. 
For  this  reason,  the  model  definition  itself  must  be  adaptive.  In  order  to  make  this  possible, 
a  hierarchic  sequence  of  models  has  been  defined.  The  essential  property  of  hierarchic  plate 
models  is  that  the  exact  solutions  corresponding  to  the  sequence  of  models  converge  to  the 
exact  solution  of  the  fuUy  three-dimensional  problem: 
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where  is  the  exact  solution  of  the  fully  three-dimensional  problem,  is  the  exact 

solution  of  the  zth  hierarchic  model.  Subsequently  the  following  desirable  properties  were 
stated  and  the  conditions  for  achieving  them  clarified: 


1. 


(a)  With  respect  to  the  thickness  {t)  approaching  zero,  the  exact  solution  of  each  model 
should  converge  to  the  exact  solution  of  the  fuUy  three-dimensional  problem  in  energy 
norm: 


lim 

t — >0 


=  0  i  =  1, 2, . . . 


In  the  case  of  plates  the  limiting  solution  is  the  Kirchhoff  plate  model.  The  first  model 
in  the  hierarchy  is  the  Reissner-Mindlin  model. 


2.  (b)  When  is  smooth  then  the  hierarchic  models  yield  optimal  rates  of  convergence. 
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where  Qj  >  is  the  largest  possible  constant. 


When  the  region  of  interest  includes  the  boundaries  of  plates  and  shells  then  boundary- 
layer  effects  must  be  taken  into  consideration.  The  exact  solution  of  plate  models  may 
differ  very  substantially  from  the  exact  solution  of  the  corresponding  fully  three-dimensional 
problem  at  the  boundaries.  Thus,  when  the  data  of  interest  depend  on  the  solution  at 
the  boundaries,  which  is  often  the  case  in  engineering  apphcations,  proper  model  selection 
through  the  use  of  hierarchic  models  is  essential. 


z/h 
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Figure  5.  Typical  through-thickness  shear  stress  distribution  computed  for  a  3-ply 
(-45  -I-  45  -  45)  soft-simply  supported  square  plate  from  a  hierarchic  sequence  of  models. 
The  reference  is  the  fully  three-dimensional  model. 

2.6  Numerical  Analysis  of  Material  Interface  Singularities  in 
Two  Dimensions 

Eigenpairs 

The  solution  for  the  hnear  elasticity  problem  in  two  dimensions  (2-D)  in  the  vicinity  of  a 
singular  point  can  be  expanded  in  the  form  [7]: 

oo  M  ^ 

W  ^  CimT^'  In’"  rflmi^)  (1) 

i=l  m=0 

where  Cim  are  the  coefficients  of  the  asymptotic  expansion  (called  the  generalized  stress 
intensity  factors  -  GSIFs),  and  ai  and  fimid)  are  eigenpairs  which  depend  on  the  boundary 
conditions. 
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At  present,  the  possibility  to  evaluate  these  GSIFs,  which  determine  failure  initiation  at 
singular  points  in  a  2-D  domain,  especially  when  the  singularity  is  caused  by  multi-material 
interfaces,  is  very  Hmited. 

A  numerical  method  based  on  the  Steklov  problem  for  the  computation  of  the  eigen- 
pairs  resulting  from  singularities  due  to  corners,  abrupt  changes  in  material  properties  and 
boundary  conditions  is  presented  (see  [8]). 

Numerical  studies  have  indicated  that  the  computed  values  converge  strongly,  are  accu¬ 
rate  and  inexpensive  from  both  the  computational  point  of  view  and  the  point  of  view  of 
human  time  needed  for  input  preparation. 

This  method  is  very  important  from  the  practical  point  of  view  because  it  provides  a 
rigorous  quantitative  basis  for  investigating  failure  events,  such  as  delaminations  of  composite 
materials  at  corners,  and  failure  in  electronic  devices. 

The  Steklov  Weak  Form. 

Notation:  We  denote  the  two  displacements  (variables)  in  the  x  and  y  directions  by  Ux 
and  Uy  respectively.  The  normal  and  tangential  displacements  and  tractions,  will  be  denoted 
by  T„,  Tt  and  tin,  Ut,  respectively.  In  the  vicinity  of  the  corner  we  assume  that  no  body 
forces  are  present. 

Let  us  consider  a  domain  CIr  shown  in  the  figure  6,  where  r,  0  are  the  coordinates  of  a 
cyhndrical  coordinate  system  located  in  the  singular  point.  On  the  boundaries  Fi  and  r2 
homogeneous  boundary  conditions  are  introduced. 


Figure  6.  Domain  and  notations  for  the  modified  Steklov  formulation. 


In  Qr,  Ux  and  Uy  may  be  represented  as  follows: 


(2) 


Under  special  (exceptional)  circumstances,  u  may  also  have  additional  terms  Inr  terms, 
however  this  case  is  not  treated  in  the  following. 

Using  (2),  on  Fs  we  have: 


{du/du)  —  {alR)u  ,  {x,y)eTz, 


(3) 
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and  a  similar  condition  on  r4. 

Multiplying  the  equilibrium  equation  by  v  =  {v^  ,  VyY  €  H^{^r)  x  we  obtain 

after  some  mathematical  manipulations  the  following  weak  form,  called  the  modified  Steklov 
weak  form: 

Seek  a  G  C  ,  0  ^  u  E  H^{Qr)  x  H^(Qr), 

B{u,  v)  +  ELi  Mi{u,  v)  -  {Mr{u,  v)  +  Ur- {u,  v))  = 
a  (Mr(u,v)  +  Mr-(u,  v)),  'iv  E  H^{^r)  x  H^{CIr)  (4) 


where: 


B{u,v)  = 

ELi  Uti(u,  v)  =  Ef=i  /r,  ^[Aif[Ak]uds, 


Mr{u,v)  =  j  ^[A\Y[Az][E][A^]ii}^^^^dd, 
e 

Ur{u,v)  =  I  [^[A^f[A,][E][D^^^^^^d9  , 


(5) 

(6) 


where  the  various  matrices  are  given  in  [9],  and  [E]  is  the  material  matrix. 

Remark  1  The  domain  Q,r  does  not  include  singular  points,  hence  no  special  refinements 
of  the  finite  element  mesh  is  required.  Furthermore,  0^  is  very  small  in  size. 

The  bilinear  forms  Ur  and  Ur-  are  non-symmetric  with  respect  to  u  and  v.  As  a 
consequence,  the  symmetric  properties  of  the  weak  form  are  destroyed.  This  means  that  in 
general  complex  eigenvalues  and  eigenvectors  exist. 

Also  note  that  the  formulation  of  the  weak  form  has  not  limited  the  domain  Qr  to 
be  isotropic,  and  in  fact  can  be  applied  to  multi-material  anisotropic  interface,  as  will  be 
demonstrated  by  a  numerical  example. 

The  expressions  in  (4)  are  reformulated  in  the  framework  of  the  p-version  of  the  finite 
element  method,  where  the  solution  of  the  eigenproblem  are  the  desired  eigenpairs. 
Numerical  Example 

Consider  two  orthotropic  materials,  graphite  and  adhesive  (epoxy),  bonded  together, 
with  plane  strain  condition  assumed.  See  figure  7. 
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Figure  7.  Orthotropic  bonded  materials  test  problem. 

This  problem  was  chosen  to  demonstrate  the  Steklov  method  for  anisotropic  multi¬ 
materials  with  a  singular  point.  The  material  properties  are  listed  in  Table  1,  where  E 
is  the  modulus  of  elasticity,  v  is  the  Poisson’s  ratio,  and  G  is  the  modulus  of  rigidity. 

The  first  three  exact  non-zero  eigenvalues,  obtained  using  the  Lekhnitskii  stress  potentials, 
are  given  in  the  following;  {ai)Ex  =  0.905  ±  O.OOOOi,  (012) ex  =  1.000  ±  O.OOOOi,  {az)EX  = 
1.944  ±  0.3051i. 

The  mesh  used  for  this  example  problem  has  the  minimum  possible  number  of  finite 
elements,  i.e.  one  element  in  each  anisotropic  subdomain.  The  exact  eigenvalues  are  given 
with  an  accuracy  of  up  to  the  third  digit,  so  that  the  accuracy  of  the  numerical  results  may 
be  assessed  up  to  about  0.01%  relative  error. 

Convergence  curves  of  the  absolute  relative  error  in  the  eigenvalues  is  shown  in  figure  8. 


OOP 

Figure  8.  Convergence  of  first  three  eigenvalues. 

The  results  demonstrate  an  excellent  convergence  rate  for  the  coarsest  mesh  possible. 


2.7  Extraction  of  GSIFs  Using  the  Principle  of  Minimum  Com¬ 
plementary  Energy 

By  utilizing  the  principle  of  minimum  complementary  energy  in  conjunction  with  the  Steklov 
formulation  and  the  p- version  of  the  finite  element  method,  the  GSIFs  for  anisotropic  as  well 
as  for  isotropic  domains,  can  be  computed  with  high  accuracy. 

The  following  method  is  proposed.  A  small  subdomain,  Qr,  is  constructed  around  the 
singular  point,  by  intersecting  a  circle  of  radius  R  centered  in  the  singular  point,  and  the 
domain  Cl.  The  approximated  pointwise  finite  element  solution  is  applied  on  the  boundaries 
of  CIr. 

The  finite  element  method  based  on  the  principle  of  minimum  complementary  energy  is 
used  now  over  CIr,  where  the  trial  and  test  functions  are  the  approximated  eigenpairs  with 
the  GSIFs  as  unknowns.  The  finite  element  formulation  chooses  these  coeflEicients  such  that 
the  complementary  energy  in  CIr  is  maximized.  Solving  the  finite  element  system  of  equation 
over  CIr,  one  obtains  an  approximation  for  the  series  coefficients. 

Numerical  experiments  for  the  scalar  elliptic  problem  (for  both  isotropic  and  anisotropic 
materials)  showed  that  the  GSIFs  converge  to  the  exact  values  as  fast  as  the  error  in  the 
strain  energy,  thus  exhibit  superconvergence.  For  a  detail  discussion  and  numerical  examples 
see  [10]. 
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Table  1:  Material  properties. 

Summary 

It  is  expected  that,  much  like  crack  extensions,  failure  initiation  and  propagation  can 
be  correlated  with  the  GSIFs,  but  neither  the  computational  procedures,  nor  the  analytical 
understanding  are  available  for  the  anisotropic  multi-material  interfaces. 

A  rehable  numerical  method  for  the  computation  of  the  eigenpairs  and  the  generalized 
stress  intensity  factors  (GSIFs)  which  characterize  the  solution  in  the  neighborhood  of  sin¬ 
gular  points  in  anisotropic  multi-material  interfaces  has  been  presented. 

Numerical  experiments  and  mathematical  analysis  indicate  that  the  computed  values 
converge  strongly,  are  accurate  and  inexpensive  from  the  points  of  view  of  human  time 
needed  for  input  data  preparation,  and  required  computer  time. 

This  method  is  very  important  because  it  provides  a  rigorous  quantitative  basis  for  inves¬ 
tigating  failure  events,  such  as  delamination  of  composite  materials,  and  failure  in  electronic 
devices. 
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